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ABSTRACT 

Full polarimetric radio interferometric calibration is performed by estimating 2 by 2 
Jones matrices representing instrumental and propagation effects. The solutions ob- 
tained in this way differ from the true solutions by a 2 by 2 unitary matrix ambiguity. 
This ambiguity is common to all stations for which a solution is obtained but it is 
different for solutions obtained at different time and frequency intervals. Therefore, 
straightforward interpolation of solutions obtained at different time and frequency in- 
tervals is not possible. In this paper, we propose to use the theory of quotient manifolds 
for obtaining correct interpolants that are immune to unitary matrix ambiguities. 

Key words: Instrumentation: interferometers; Methods: numerical; Techniques: in- 
terferometric 



1 INTRODUCTION 

Most modern radio interferometers have dual polarized feeds 
and therefore, the us e of the matrix measurement equation 
(|Hamaker et al.|[l996l ) gives a compact and accurate descrip- 
tion of their operation. Calibration of such an interferometer 
essentially boils down to the estimation of Jones matrices 
of size 2x2 with complex entries. As shown by lHamakeil 
(|2000i ) . the solutions acquired for the Jones matrices will 
only be equivalent to the true solution upto a unitary ma- 
trix ambiguity. 

This ambiguity would not hinder further processing of 
data because it cancels out during correction of the data 
using the obtained solutions. However, the ambiguities do 
prevent us from using the solutions for further modelin g 
of instrumental effects (e.g., beam shape (|Yatawattall2012l )) 
and eff ects due to the prop agation medium such as the iono- 
sphere (|lntema et al.ll2009l ). Moreover, interferometers oper- 
ating at low frequencies have a wide field of view and cali- 
bration has to be pe rformed along hund reds of directions in 
an efficient manner (|Kazemi et al .1120111 '). Along each direc- 
tion, we would have an ambiguity which is independent of 
other directions. 

In this paper, we present a method of interpolation of 
the calibration solutions (Jones matrices) and we consider 
the case where each solution is affected by any unknown 
unitary ambiguity. We consider the simplest case of interpo- 
lation to present our method: Finding the mean of a given 
set of solutions which can also be extended to interpolation 
with weighted averaging. The uses of interpolation or aver- 
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aging are numerous. First, averaging of solutions obtained at 
adjacent time and frequency intervals provides us with a ro- 
bust estimate of calibration solutions especially under noisy 
situations. Moreover, interpolation reduces the number of 
data points that needs to be visualized. This is important 
when we have solutions over hundreds of directions in the 
sky at a number of time and frequency intervals. Interpo- 
lation also provides us solutions when there is missing or 
flagged data points. 

Traditional calibration software such as AIPS that are 
based on a scalar data model has the ability to interpolate 
scalars such as the amplitude or phase of a single polar- 
ization. However, this is not possible when we are dealing 
with 2x2 matrices as what is presented in this paper. Due 
to the unitary ambiguity, any linear operation in Euclidean 
space with such matrices would not give us a feasible inter- 
p pla nt. Therefore, we explore the quotient manifold struc- 
ture (|Absil et al. I [200^ 1 of the calibration solutions. We use 
the theory of interpolation over manifolds to get a feasible 
solution to our problem. We present an algorithm to find 
the mean of a given set of calibration solutions. This al- 
gorithm could also be extended to interpolation with any 
positive weighting scheme. The interpolant we obtain using 
this algorithm still has a unitary ambiguity but it is closer to 
the true interpolant (without ambiguities) than what we ob- 
tain by standard (Euclidean) interpolation. Similar work on 
averaging or interpolation over manifolds and Lie groups ap- 
pear in diverse a reas o f resea r ch and we refer the reader t o 
iFiori fc Tanakal J2009h. [Fioril (l201ll ). iKaneko et all (|2012h . 
and I Amsallem fc Farhat |2008l ) for further information. 

The rest of the paper is organized as follows: In section 
[2] we give an overview of radio interferometric calibration 



© 2012 RAS 



2 Sarod Yatawatta 



and the ensuing ambiguities. In section [3] we present the in- 
terpolation scheme. We provide simulation results to prove 
the superiority of the proposed scheme in section [4] and fi- 
nally, draw our conclusions in section [5] 

Notation: Matrices and vectors are denoted by bold up- 
per and lower case letters as J and v, respectively. The trans- 
pose and the Hermitian transpose are given by(.) T and(.) H , 
respectively. The matrix Frobenius norm is given by ||.||. The 
set of complex numbers is denoted by C. The identity ma- 
trix is given by I. The angle of a complex number is given 
by Z. 



2 RADIO INTERFEROMETRIC 
CALIBRATION 

Consider a radio interferometer with N stations. Let the 
sky signal consist of M discrete sources. The observed data 
on the baseline pq fo rmed by stations p and q is given by 
jHamaker et alJll996h 

M 

+ N P9 (1) 

m— 1 

where V P9 and N pq (€ C 2x2 ) are the visibility matrix, and 
the noise matrix, respectively. The instrumental and propa- 
gation parameters of stations p and q along the direction of 
the m-th source are represented by J pm and 3 qrn (6 C 2x2 ), 
respectively. The source coherency matrix of the m-th source 
is given by C pgm (G C 2x2 ). Throughout this paper, we as- 
sume all the sources are unpolarized and therefore, C pqm 
to be diagonal. Calibration is the esti mation of J pm fo r all 
p G [I, N] and m G [1, M]. As noted bv lHamakerl (|2000l ). for 
any unitary matrix U m (G C 2x2 ), U*U m = U m U„ = I, a 
valid calibration solution would also be J pm U m . Note that 
this unitary ambiguity U m can have different values for dif- 
ferent 77i, or different directions. 

Even with a solution that has an ambiguity, the data 
could still be corrected for the estimated errors because 
the ambiguity will cancel out during correction. However, 
the ambiguity prevents us from using the calibration solu- 
tions to model the instrument (such as the beam shape) 
or propagation phenomena (such as the ones that happen 
in the ionosphere). The simplest example of further use of 
solutions is finding their mean. Let us consider having two 
solutions for station p along the direction m, say at adja- 
cent time or frequency intervals: J pm iU m i and J pm 2U m 2. 
Also assume that the true values of the Jones matrices 
are almost identical, i.e., J pm i ~ J pm 2 = J P m- The true 
mean is (J pm i + J pm 2)/2 = J pm while we obtain the in- 
terpolant (J pm iU m i + J pm 2U m2 )/2 = J pm (U m i +U m2 )/2. 
Due to the fact that (U m i + U m 2)/2 is not a unitary matrix 
(U^ 1 U m 2 7^ 0), the interpolant is not a valid solution for 
J pm anymore and will not satisfy {TJ. 

We can rewrite |T]) as 

M 

^P<3 = / . A p J m C pgm J m A^ + N P9 (2) 

m— 1 

where J m is the augmented matrix of all Jones matrices 
along the m-th direction 

Jra= [Jlmr-'jJwm] , GC . (3) 



The canonical selection matrix A p (and A q likewise) is given 
as 

A p = [0,0,...,I,...,0], GC 2x2iv . (4) 

In Q, the p-th block of columns is a 2 x 2 identity matrix 
while the rest is all 0. Then, we see that J m U m (where U m 
is unitary) is a valid solution for J m . 

In this paper, we try to solve the following problem: 
Given a set of solutions whose intrinsic values (i.e., solutions 
without an ambiguity) are almost equal, we find the inter- 
polant that is immune to unitary ambiguities. Let the set 
of solutions (taken at adjacent time and frequency intervals 
and even along adjacent directions) be 5, 

5 = {Ji,J 2 ,...,Jk} (5) 

that has K elements. The elements in S satisfy 

J fe = J fc U fe , k G [l,K] (6) 

where is the intrinsic value of solution and is the 
unitary ambiguity. We make the additional assumption that 
all intrinsic values are almost the same, i.e., 

Ji^J 2 ...~Jk. (7) 

This assumption mostly holds for solutions obtained at ad- 
jacent time and frequency intervals as well as along adjacent 
directions, provided that the scale difference due to source 
models along adjacent directions is taken into account. We 
make this assumption so that we are not affected by any 
aliasing errors. Note that the unitary matrices are almost 
surely not equal 

Ui^U 2 ...^U K . (8) 

We restate the problem we need to solve: Given the set 
of solutions 5, find the mean J of the solutions such that it 
is the most accurate approximation of the intrinsic mean 

J=~£J* 0) 

fe=i 

within a unitary ambiguity. By replacing the summation 
with weighted summation in ([9]), this could also be extended 
to interpolation with any positive weighting scheme. 



3 INTERPOLATION 

In this section, we first present the well established scalar 
averaging technique used in radio interferometry. As men- 
tioned before, this does not extend to the case with Jones 
matrices as calibration solutions. In order to proceed further, 
we give a rather pedagogical overview of the manifold geom- 
etry of calibration solutions. Next, we present an averaging 
algorithm (Algorithm I) using the quotient manifold in sec- 
tion 13.31 We also give an alternative algorithm (mainly for 
comparison purposes) that averages using the tangent space 
of the quotient manifold (Algorithm II) in section 13.41 
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3.1 Interpolation: The scalar case 

The calibration solutions of N stations for an interferometer 
with a single polarization can be given as 



g = 



91 




92 








. 9N _ 





l<7i|e J 



\9N\e 



(10) 



where ip is the phase ambiguity common to all stations. Con- 
sider the averaging of K such solutions given by the set Q 
as 



Q = {gie JVl ,g2e- 



(11) 



Consider the calculation of the average phase for station 
n using the solutions in the set g. Normally, we keep one 
station (out of N) as the reference (say the first station). 
With one station kept as the reference, the average phase of 
station n becomes 

i K 

Zg„ = — ^2{Zg nk + ipk - ^gik - iph) (12) 

k=l 

1 K 1 K 

k=l k=l 

where J^fcLi ^-9nk is the intrinsic average. Moreover, the 
term i Y^=i z 9 lk i s common to the averaged phases of all 
stations. Therefore, phase averaging can be solved for the 
single polarized case upto a common phase ambiguity. As 
the ambiguity e 3 * does not affect the amplitudes of the so- 
lutions, amplitudes can also be averaged without hindrance. 
Furthermore, by using positive weights in the summation of 
(|12p . the same method can be applied to any interpolation 
scheme. This form of averaging and interpolation is widely 
used in current interferometric data processing. However, 
this method does not extend to calibration solutions with 
dual polarized interferometers where we have Jones matri- 
ces as our solutions. 



3.2 The quotient manifold structure of calibration 
solutions 

We provide a brief overview of manifolds and Lie groups 
before we proce ed. A mor e gen eral overview of th is subject 
can be found in (|201lh and lAbsil et~ai1 (|2008h . A man- 
ifold can be loosely described as a set of entities, together 
with a set of mappings (charts) that can locally describe the 
manifold in Euclidean space. A "quotient" manifold is a sub- 
manifold of a larger manifold and the entities in the quotient 
manifold represent more than one entity in the embedding 
manifold. 

This notion of a quotient manifold naturally represents 
the calibration solutions with unitary ambiguities. Given the 
set S in (J5]), we consider two solutions (say Ji and J2) to be 
"similar" if they are related by a unitary matrix, i.e., 

Ji ~ J 2 ^ Ji = J2U12 (13) 

where U12 is unitary. The equivalence relation ~ satisfies 
re flexive, symmetric , and transitive conditions as described 
111 lAbsil et al.l (|2008f ) . Therefore, assuming the intrinsic value 
of all elements in S are the same (|7J, we can select only 




Figure 1. The quotient manifold geometry l|Absil et alj|2008h of 
the calibration solutions. The dashed (blue) line (on A4) repre- 
sents the equivalence class of all solutions that are related to J 
by a unitary ambiguity. This equivalence class is represented by 
a single point on the quotient manifold M — M.J ~. The vertical 
space Vj is the vector space tangential to the equivalence class 
and the horizontal space Hj is the orthogonal complement. 

one element from S to represent the whole set, under the 
equivalence relation ~, given by (|13[) . 

Consider Ai to be the manifold of all 2N x 2 complex 
matrices (C 2JVx2 ), then, we can represent all elements in S 
on the quotient manifold M = M/ ~, where the equiva- 
lence relation is given by (|13jl. The mapping n (canonical 
projection) is defined such that any matrix JU on M (U 
unitary) is mapped onto a single point, 7r(J) on M/ ~. 

With the mapping tt, we have the equivalence class 

7r _1 (7r(J)) = {JU : UU fl = U H U = I,U G C 2x2 } (14) 

of solutions that is represented by a single point on M/ ~, 
as illustrated in Fig. [T] 

The vertical space and horizontal space are vector 
spaces that are related to the manifold as follows. We take 
the vertical space to be the tangent space to the equivalence 
class 7r _1 (7r(J)) at J, i.e., 

. h 



Vj = {JA : A' 



-A, A € C 2X2 } 



(15) 



and we choose the orthogonal complement of the vertical 
space as the horizontal space %j, 



Hj = {J±B : B <= 



V-2X2 j 
~i2Nx2N- 



(16) 

The orthogonal complement (G C 2iv X2iv ~ 2 ) is a matrix 
whose columns are orthogonal to those of the matrix J, i.e., 
J H J ± = 0. 

With this formal representation of the manifold geom- 
etry of the solution space at hand, we present two algo- 
rithms for interpolation in sections l3.3l and 1 3.41 Similar aver- 
aging techniques on o ther manifolds are alr eady being inves- 
tigated. For instance, iKaneko et all (|2012T ) present an algo- 
rithm for averaging on a compact Stiefel manifold. A similar 
al gorithm for averaging on th e Grassmann manifold is given 
by Amsallem fc Farhatl (|2008l ) . While a manifold has a more 



geometric structure, a Lie group has a more algebraic struc- 
ture. A Lie group can be described as a set of entities, with 
an identity element and operations for multiplication and 
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inverse (ffu]|201ll ). There is a close relation between smooth 
manifolds and Lie groups and in lFiori fc Tanakal |2009l ). sev- 
eral averaging techniques for sq uare matrice s using the Lie 
group structure is presented. In iFioril | 20 111 ), the Lie group 
structure as well as the manifold geometry is exploited for 
averaging. 

The algorithm presented in section 13.31 (Algorithm I), 
performs the averaging directly on the quotient manifold 
while the algorithm presented in section [3~4l (Algorithm II), 
performs the averaging in the tangent space to the mani- 
fold. Algorithm II is presented mainly for comparison with 
Algorithm I, as mo s t existing work on averag ing, such as 
(|Kaneko et al.ll2012l ; lAmsallem fc Farhatl 120081 ), is done in 
the tangent space. Both these algorithms can be extended 
to interpolation by using any positive weighting scheme. 



3.3 Averaging on the quotient manifold 
(Algorithm I) 

We present the algorithm to find the mean of the set iS in 
([5j|. Let us call the estimated mean as J. The basic idea is 
to find the element on the quotient manifold that represents 
the set S as accurately as possible. This would also be the 
mean of the set S. 

(i) Start with the initial estimate as one value from S, 
say J <— Ji. The error threshold is given by e. 

(ii) For each element in S, find unitary Pfc (£ C 2x2 ) such 
that 



arg mm 
p«p i .=p l p«=i 



JfePi 



(17) 



This is basically an "alignment" operation and the details 
of this step are given in section [331 
(iii) Form 



1 K 



(18) 



(iv) Find the unitary projector P to minimize || J — GP|| 2 
as given in section 1531 

(v) If || J - GP|| < e then stop. Else update J <- G and 
go to step (ii). 

(vi) Return J as the mean. 

Note that this algorithm can be modified for interpo- 
lation by replacing the averages in (|18|) with weighted av- 
eraging using positive weights. The proof of convergence of 
this algorithm is part of future research. For the moment, 
we rely on numerical simulations to test its convergence in 
section [4] 



3.4 Averaging in the tangent space (Algorithm II) 

The algorithm presented in this section projects each ele- 
ment in the set S to the horizontal space of the quotient 
manifold before averaging is performed. As before, let us 
call the estimated mean as J. 

(i) Start with the initial estimate as one value from S, 
say J <— Ji. The error threshold is given by e. For better 
convergence, a positive scalar p € (0, 1] is used. 



(ii) Form the orthogonal projector matrix V (£ (£ 2iVx2JV ) 

V = I-1(J H J) 1 J H . (19) 

(iii) Project each element in <S, onto the horizontal space 



T~L~t as 



W fc = VJ fe . 
(iv) Form the average in Hj as 



(20) 



(21) 



(v) Form current estimate for the average as G = J+pW. 
Find the unitary projector P to minimize || J— GP|| 2 as given 
in section [531 

(vi) If ||J — GP|| < e then stop. Else update J <— G and 
go to step (ii). 

(vii) Return J as the mean. 

3.5 Finding unitary projector to minimize 

II J- JfcPfcll 2 

Wha t we have to solve i s in fact the matrix Procrustes prob- 
lem (ISchoneman"nlll966l ) and we use the algorithm given in 
IHighaml (|200Sf ). 



(i) Find the product X = j" J. 

(ii) Find the singular value decomposition of X as 

USV H = X. (22) 

(iii) Return P fc = UV H . 

The proof can be found in IHighaml (|2008l ) . 



3.6 Discussion 

The main assumption used in Algorithm I is that all ele- 
ments in S belong to the equivalence class (dashed line in 
Fig. [1} or are very close to it. So on the quotient manifold, 
they lie on a small area that can be considered locally Eu- 
clidean, thus enabling averaging. For Algorithm II, we per- 
form the averaging in the tangent space, which is a vector 
(Euclidean) space and therefore averaging works. Computa- 
tionally, Algorithm II is much more expensive because of the 
calculation of the projection matrix at each iteration. Also, 
it is numerically less stable and hence the need of p for bet- 
ter convergence. We use Algorithm II mainly for comparing 
the performance of Algorithm I. The fundamental question 
posed here is whether (for our specific problem) it is bet- 
ter to average on the quotient manifold or in the tangent 
space. Most existing work on averaging or interpolation use 
operations in the tangent space therefore comparison of Al- 
gorithms I and II in terms of their performance is important. 
In the next section, we give simulation results to compare 
their performance. 



4 SIMULATION RESULTS 

We consider an interferometric observation with N = 40 
stations and interpolation of K = 10 solutions. Therefore 
the set S in ((5]) has cardinality of 10, with each matrix Jfc 
of size 80 x 2. The intrinsic values of the first matrix, Ji, 
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are generated as U(0, 1) + jU(0, 1) (drawn from a uniform 
distribution in [0, 1]). The intrinsic values of the remaining 
matrices J2, . . . , Jk are obtained by perturbing the elements 
of Ji by adding a(U(0, l)+jU{0, 1)) to them. The value for 
a is varied for different simulations as explained later. How- 
ever, we keep the value of a in [0, 0.5] to ensure the intrinsic 
variation is small enough that averaging (or interpolation) 
is not dominated by aliasing. Once the matrices are gener- 
ated in this fashion, we multiply each intrinsic matrix Jfe by 
a random unitary matrix Ufe to get Jfc = JfcUfc. To each 
realization of Jfe, a noise matrix N is added. The elements 
of N are generated to be complex circular Gaussian random 
variables. The values of N are scaled to get a given signal 
to noise ratio (SNR) where the SNR is defined by 



SNR: 



IINI 



(23) 



before adding them to Jfe. 

We also calculate the normal (or Euclidean) average for 
comparison, as 



1 K 



K 



(24) 



Moreover, we also calculate the intrinsic sample variance, 
using the intrinsic mean of @ as 



var(J) 



K||J|| 2 El 



Eh j - j kI 



(25) 



The criterion that we use for measuring the performance 
of the various averaging algorithms is the normalized error 
(NE), defined as 



NE 



A || J -_JU|| 
IIJF 



(26) 



where J is the intrinsic mean, J is the estimated mean, using 
(i) Algorithm I (ii) Algorithm II or, (iii) Euclidean average 
of (|24[) and finally, U is a unitary projector determined as 

in section [3~5l to minimize IIJ — JUII 2 . 



Normal 
Algorithm I 
Algorithm II 
Sample Variance " 



g 0.03 - 




40 50 60 
No. of realization 



Figure 2. Normalized error for 100 realizations of <S with <r = 
0.1 and SNR = 100. The intrinsic sample variance is given by 
the dashed line. The proposed Algorithm I performs much better 
than normal (Euclidean) averaging and almost at the level of 
intrinsic variance. Algorithm II performs better than Euclidean 
averaging but is worse than Algorithm I. 




4.1 Simulation I 

We generate the set S as described above 100 times, keeping 
a — 0.1 and SNR = 100. For each realization, we estimate 
the average by Algorithm I, using 10 iterations and Algo- 
rithm II, using 60 iterations (the reason for this will be ex- 
plained later) with p = 0.3. We have shown the normalized 
estimation error (NE) in Fig. [2] for different approaches. We 
have also plotted the intrinsic variance, calculated using (|25p 
in this figure. As seen in this figure, the proposed Algorithm 
I has an error almost comparable with the intrinsic variance. 
The proposed Algorithm II has a slightly worse performance 
but it is still better than Euclidean averaging. 

In Fig.[3]and Fig. [4] we have shown the convergence per- 
formance of both Algorithm I and Algorithm II. We measure 
the convergence by the norm of the difference between the 
current estimate and the updated estimate at each iteration. 
It is clear that Algorithm I has much better convergence 
(only about 4 iterations) than Algorithm II, which does not 
converge even after 60 iterations. 



Figure 3. The convergence of Algorithm I for 100 different real- 
izations. After about 4 iterations, no further improvement occurs 
and therefore, we can fix the maximum number of iterations at 
4. 



4.2 Simulation II 

We vary the values of a and SNR in this simulation. For 
each value of a and SNR, we generate 100 realizations of 
the set S as before and compute the average. In Fig. [5l 
we have shown the mean normalized error over all realiza- 
tions for both proposed algorithms as well as for Euclidean 
averaging. For Algorithm I, we used 4 iterations and for Al- 
gorithm II, we used 60 iterations with p = 0.3. From Fig. [S] 
it is clear that Algorithm I gives the best results. Moreover, 
the dependence on a is not significant. The performance of 
all three methods degrade at low values of SNR. Algorithm 
II performs worse than Algorithm I and comparing the com- 
putational cost, Algorithm I clearly gives better results. 
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Figure 4. The convergence of Algorithm II for 100 different 
realizations. Even after 60 iterations there are some cases that 
do not converge. 
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Figure 5. The mean normalized error over 100 realizations for 
various values of cr and SNR. Algorithm I performs significantly 
better than normal (Euclidean) averaging while Algorithm II 
performs slightly better but worse than Algorithm I. Algorithm 

I is almost insensitive to the values of cr, especially at high values 
of SNR. 

4.3 Simulation III 

In this simulation, we perform weighted averaging, where 
the weights are generated from a uniform distribution in 
[0, 1]. Once again, for each value of a and SNR, we generate 
100 realizations of S and for each realization, we generate 
a new set of weights. We emphasize that Algorithm II did 
not give convergent results even at 60 iterations and for any 
value of p. Therefore, we omitted the results of Algorithm 

II because it clearly fails to perform weighted averaging. 



5 CONCLUSIONS 

We have presented a method for averaging Jones matri- 
ces obtained in radio interferometric calibration by exploit- 



0.14 
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Figure 6. The mean normalized error over 100 realizations for 
various values of a and SNR under weighted averaging. Algo- 
rithm I performs better than Euclidean averaging while Algo- 
rithm II fails to perform well. 

ing the quotient manifold structure of such solutions. This 
method could also be used for weighted averaging, or in- 
terpolation. Unlike Euclidean averaging, which gives inac- 
curate results due to the unknown unitary ambiguity in the 
solutions, the proposed methods give significantly better re- 
sults. For comparison, we have also proposed an alternative 
algorithm that operates in the tangent space to the mani- 
fold. Simulation results suggest that averaging directly on 
the quotient manifold is better than averaging by using the 
tangent space, for our specific problem. One possible reason 
for this could be due to numerical instability. Existing work 
that exploit the tangent space for averaging have matrices 
with orthonormal columns (e.g. the Stiefel manifold). How- 
ever, in our problem we do not have such a constraint and 
this could result in numerical instability. Future work would 
focus on improving the numerical stability and interpolating 
along geodesies of the quotient manifold. 
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